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Abstract. We present a rigorous scheme that makes it possible to compute eigen- 
values of the Laplace operator on hyperbolic surfaces within a given precision. The 
method is based on an adaptation of the method of particular solutions to the case 
of locally symmetric spaces and on explicit estimates for the approximation of cigen- 
functions on hyperbolic surfaces by certain basis functions. It can be applied to check 
whether or not there is an eigenvalue in an e-neighborhood of a given number A > 0. 
This makes it possible to find all the eigenvalues in a specified interval, up to a given 
precision with rigorous error estimates. The method converges exponentially fast with 
the number of basis functions used. Combining the knowledge of the eigenvalues with 
the Selberg trace formula we are able to compute values and derivatives of the spec- 
tral zeta function again with error bounds. As an example we calculate the spectral 
determinant and the Casimir energy of the Bolza surface and other surfaces. 

I. Introduction 

Let (M, g) be a compact Riemannian manifold of dimension n and let A be the 
(positive) Laplace operator on functions, which in local coordinates is given by 



Here \g\ denotes the determinant of the metric tensor and g are the components of 
the dual metric on the cotangent bundle. This operator is formally self-adjoint on the 
space of smooth functions with inner product 



defined by the Riemannian measure \J\g\dxi ■ ■ ■ dx n . Then A extends to a self-adjoint 
operator L 2 (M) D H 2 (M) — > L 2 (M) with compact resolvent. This means there exists 
an orthonormal basis {4>j \ j G N } in L 2 (M) consisting of eigenfunctions 





dxi, ■ ■ ■ dx, 



n 




A0j = \j(f>j 



which we assume to be ordered such that = A < Ai < A2 < . . .. 
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It is a classical problem to compute the eigenvalues and the eigenfunctions of the 
Laplace operator on a given manifold. Its solution allows complete control over the 
functional calculus and thus over solutions to the heat equation, the wave equation, or 
the Schrodinger equation. 

In this paper we show how an adaptation of the method of particular solutions can 
be used to compute eigenvalues on manifolds with a high accuracy, beyond of what can 
currently be achieved by finite element methods or boundary element methods. We will 
focus primarily on the case of two dimensional oriented surfaces of constant curvature. 
These are the simplest examples of manifolds with non-trivial spectral geometry. They 
are topologically classified by their genus. By the Gauss-Bonnet theorem the curvature 
k of a surface of genus g can always be normalized to be k = sign(g — 1) by multiplying 
the metric by a positive constant. This divides oriented surfaces of constant curvature 
into three categories 

• k — +l,g = 0: up to isometries there is only one such manifold: the round 
sphere S 2 . The spectrum of the Laplace operator on the round sphere is of 
course well known and the spherical harmonics provide a basis consisting of 
eigenfunctions. 

• k — 0,g = 1: such manifolds are isometric to flat tori which are obtained as 
quotients of IR 2 by co-compact lattices. The moduli space of equivalence classes 
of flat metrics on a given topological torus is the modular surface SL(2,Z)\H, 
where EI = {x + iy G C | y > 0} is the upper half space. For a surface M?/L 
a basis consisting of eigenfunctions is obtained by taking the Fourier modes 
associated with points in the dual lattice L* and the corresponding eigenvalues 
are the squares of the norms of these points. 

• k — —1, g > 1: these are the so called hyperbolic surfaces. They can be obtained 
either as quotients r\H by co-compact hyperbolic lattices in SL(2,R) or by 
glueing hyperbolic pairs of pants. The moduli space of hyperbolic metrics on a 
given topological surface is a quotient of the Teichmiiller space T g by a discrete 
group, the mapping class group. The Teichmiiller space for a genus g surface has 
dimension 6g — 6 so that a constant curvature metric on a topological surface 
may be given by specifying 6g — 6 parameters. Unlike in the previous two cases 
the spectrum of the Laplace operator of a hyperbolic surface can not currently 
be explicitly computed. 

The main purpose of this article is to establish an algorithm to compute the eigenval- 
ues and eigenfunctions of the Laplace operator on a hyperbolic surface up to a certain 
precision and with mathematically rigorous error bounds. Once this is achieved for a 
number of eigenvalues it is then possible to compute values of the spectral zeta function 
and the spectral determinant of the Laplace operator on such surfaces. We demonstrate 
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how this can be done using the Selberg trace formula in such a way that the error can 
again be bounded rigorously. 

The method we use is adapted to a pants decomposition of the surface. For a par- 
ticular basis of functions we establish various bounds and estimates that prove that 
the method of particular solutions can be applied to an interval and yields all eigen- 
values in that interval. In particular we prove a novel estimate that allows us to show 
non-existence of eigenvalues in a certain interval. 

1.1. Organization of the article and results. In section [2] we describe how an adap- 
tation of the method of particular solutions can be applied to manifolds. The method 
uses a finite dimensional space of test functions that solve the eigenvalue equation with 
eigenvalue A on open submanifolds that are glued together along co-dimension one 
hypersurfaces. The glueing conditions for the eigenfunctions, continuity of the eigen- 
functions and their normal derivatives, is measured in terms of the differences of the 
function values and normal derivatives along the hypersurfaces. For a test function 
we introduce in Equation ([T| the number F_ 1/2,-3/2(0) that measures these differences 
in suitable Sobolev norms. We prove that if ||0||l 2 (m) = 1 and -F-1/2, -3/2(0) is small 



then A must be close to an eigenvalue. Our Theorem 2.1 together with the subsequent 
estimates on the constants gives a quantitative version of this statement for hyperbolic 
surfaces that are glued along geodesic segments. 

In section [3] we review the construction of hyperbolic surfaces from pairs of pants. 
This gives an explicit parametrization of Teichmiiller space in terms of Fenchel-Nielsen 
coordinates and shows that hyperbolic surfaces can be cut open along geodesic segments 
into pairs of pants and subsequently into subsets of hyperbolic cylinders. 

In section [4] we describe a set of basis functions for a surface of genus g that is 
decomposed into pairs of pants. We prove that true eigenfunctions can be approximated 
by linear combinations of this set of basis functions with explicit bounds on the error. 



Section 4.3 deals with the construction of m x k matrices B°, Ba and Ca with the 
following properties. 

(1) The distance of A to the spectrum can be bounded from above in terms of the 



first singular value o"i(B°) of B° and its singular vector (Theorem 4.5). 



(2) The smallest relative singular value cti(Ba, Ca) is bounded from above by 

ci(k, A) + c 2 (A)dist(spec(A), A), 
where Ci(k, A) and 02(A) are explicitly computable constants and c\ is exponen- 



tially decaying in k (Theorem 4.4). 



Whereas the first property allows to prove that an eigenvalue is in a certain interval, 
the second property allows to determine intervals in which there are no eigenvalues. 
Both estimates together can be used to find all eigenvalues in a specified interval. We 
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demonstrate that these matrices can be computed within a given precision and show that 
the singular values can be bounded from above and below using interval arithmetics. As 
a proof of concept we implemented our method in Fortran and in Mathematica. This 
resulted in programs that allow to compute eigenvalues rather accurately for a surface 
of genus g with given Fenchel-Nielsen coordinates. A Mathematica program was used 
to compute the first eigenvalues of the Bolza surface with extremely high accuracy. 

The spectral zeta function (a{s) is defined as the meromorphic continuation of the 
function 

oo 
i=l 

where (Aj), 6 N are the non-zero eigenvalues of the Laplace operator repeated according 
to their multiplicities. The (zeta- regularized) spectral determinant det^(A) is defined 
by logdet^(A) = — Ca(0) (zero is not a pole of the meromorphic continuation). Since 
the meromorphic continuation is contructed from the full spectrum it is a priori not 
enough to know a finite part of the spectrum to compute the spectral determinant up 
to a certain accuracy. However, we show in section [6] that the Selberg trace formula 
may be used in conjunction with the list of eigenvalues up to a certain threshold to 
calculate values of the spectral zeta function, in particular the Casimir energy and 
the spectral determinant, up to a certain precision depending on that threshold. The 
spectral determinant for the Laplace operator on hyperbolic surfaces is of particular 
importance. In the seminal paper |OPS88] Osgood, Phillips and Sarnak showed that 
the spectral determinant is maximized at the hyperbolic metric in each conformal class 
of a given volume. Thus, understanding the extremal properties of the determinant as 
a function on the space of metrics of fixed volume in dimension two is equivalent to 
the understanding of the spectral determinant as a function on the Teichmiiller space 
of hyperbolic metrics. Of course values of the determinant for non-hyperbolic surfaces 
can be computed from the value for the corresponding uniformized hyperbolic surface 
using the Polyakov formula. 

Finally, section [7] contains of collection of examples of interesting surfaces of genus 
two and three for which we computed the spectral determinant and the value of the 
spectral zeta function at the point —1/2. The set of examples contains the isolated 
surfaces of genus two with large symmetry group which are well known to be critical 
points for any value of the spectral zeta function or the spectral determinant (see e.g. 
[KKK09J). We conjecture that the value of the spectral determinant is maximized at the 
Bolza surface and we compute the value of the spectral determinant rather accurately. 
Such a conjecture would imply that the global maximum of the spectral determinant 
as a function on the space of all metrics of fixed volume is attained at this point. The 
appendix contains explicit estimate for various resolvents as well as explicit estimates of 
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the L°° norm of eigenfunctions and derivatives of eigenfunctions. These estimates are 
needed to make the constants in our estimates explicit but they may also be interesting 
in their own right. 

1.2. Discussion. Eigenvalues of hyperbolic surfaces of genus two have been calculated 
in the physics literature by Aurich and Steiner in |AS89j using the finite element method, 
and in |AS93j using the boundary element method. In both cases the authors relied 
on the realization of the surface by geodesic octagons. It is quite interesting that the 
Hadamard-Gutzwiller model discussed in |AS89j is actually the same as the Bolza sur- 
face. Using the group action and the decomposition into geodesies triangles Ninnemann 
|Nin95j computed the eigenvalues of the regular geodesic octagon with every second 
side identified. The resulting surface, which he also refers to as the Gutzwiller octagon, 
does however not coincide with the Bolza surface. Its Fenchel Nielsen coordinates are 



given in Section 7.4 There as well as in other parts of the physics literature the Bolza 
surface is referred to as the regular octagon. 

The method of particular solutions is based on an article by Fox, Henrici and Moler 



|FHM67j . It was subsequently further developed and revived by Betcke and Trefethen 
[BT05] to achieve high accuracy in eigenvalue computations on domains with Dirichlet 
boundary conditions. Further versions, including domain decompositions and matching 
of boundary data including the normal derivative were described in IBB TO] and |Bet07j . 
Explicit error bounds for the eigenvalues in the method of particular solutions for do- 
mains in W 1 were given in |MP68j and further improved in by Alex Barnett |Bar09] . We 
would also like to refer to the latter article for further background and literature on the 
method of particular solutions. Our estimate for the error bound does not use jMP68j 
but is instead based on the more refined information contained in the resolvent of the 
Laplace operator on the manifold. Thus, for generic eigenfunctions it is expected to 
give an improvement of the order of the square root of the eigenvalues. Estimates that 
provide such an improvement without this genericity assumption have recently been 
obtained for domains in R n in |Bar09j and |BHllj . 

Finally, for hyperbolic surfaces with cusps Hejhal ( |Hej92| and |Hej99| ) introduced 
a method to compute Maass cusp forms which is based on the group action and the 
expansion of the embedded eigenvalues on the cusp. After most of the work in this 
paper was completed we learned that Booker, Strombergsson and Venkatesh [BSV06J 
have recently used the pre-trace formula together with the Taylor expansion of boundary 
data of quasi-modes to rigorously verify embedded eigenvalues for the modular surface 
SL(2, K)\H. Their way of certifying eigenvalues is in spirit similar to the method we use 
to prove that our computed eigenvalues are accurate within the error bounds. Apart 
from the fact that their method applies to a different geometric situation (surfaces 
with at least one cusp) the method of bounding the error is different when it comes to 
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technical detail: for example in |BSV06j it is estimated in terms of the L°°-norm of the 
boundary data instead of the L 2 -norm. It would be very interesting in the future to 
combine these ideas. In particular a modification of our estimates applies to hyperbolic 
surfaces with cusps. The problem of adjusting the step size in the search for Maass cusp 
forms on SL(2, Z)\H in such a way that no eigenvalues are missed could be tackled in 
this way. 

Formulae for the spectral determinant based on the length spectrum were given in 
the mathematics literature by Fried [Fri86j and by Pollicott and Rocha [PR97J, and in 
the physics literature by Aurich and Steiner [St87L IAS92] . In |PR97] , exponential con- 
vergence was proved and a numerical algorithm was established for the case of surfaces 
of genus 2 in mw-Fenchel Nielsen coordinates without twisting. We are not able to 
confirm the numerical values obtained in |PR97j in the three examples there but obtain 
quite different values. We believe that the part of the length spectrum computed in 
|PR97j was not sufficient to obtain the correct values. In contrast to the other methods 
employed our approach to compute the spectral determinant and the values for the zeta 
function allows for explicit error estimates even if the length spectrum is unknown. 

2. The method of particular solutions on manifolds 

The method of particular solutions is a method to approximate eigenvalues and eigen- 
functions of a differential operator. We will consider here the case of the Laplace op- 
erator A : C°°(M) — > C°°(M) acting on functions on a compact oriented Riemannian 
manifold M. 

To fix notations suppose that r C M is a closed subset which is the finite union 
of oriented compact codimension one submanifolds r a , possibly with boundary. We 
will assume here for simplicity that the T a intersect at most at boundary points, i.e. 
T a r\T/3 C dT a (IdTp. Removing T from M results in an open manifold M\T that is the 
interior of a manifold with piecewise smooth boundary. This manifold might however 
have corners or other singularities. We would like to define the class of functions that 
are "smooth up to the boundary" on M\T. In order to do this let us be more precise 
about how the manifold with piecewise smooth boundary is constructed from M and T. 
The space M\T equipped with the geodesic distance is not complete and we denote by 
M\r its abstract metric completion. Since M is complete there is a natural continuous 
surjection 7r : M\T — > M. We define a function / e C(M\T) to be smooth if and only 
if for every point x e M\F there exists an open neighborhood U such that there exists 
a smooth function g on M with f\u = n*{g)\u- The corresponding sheaf of smooth 
functions endows M\T with the structure of a smooth manifold with piecewise smooth 
boundary. We write C°°(M\T) for the space of smooth functions on M\T. 



Of course an element <fi G C°°(M\r) can be understood as a function in L°°(M) as 
the boundary of M\T has zero measure and the interior of M\T can be identified with 
M\T. We will now make this identification without further mention. Since each T a 
is oriented there is a natural unit normal vector field n a to T a and a closed tubular 
neighborhood of r Q diffeomorphic to T a x [— e, e]. A function G C°°(M\r) therefore 
has two boundary values on T a , a right boundary value 0+ and a left boundary value 0~. 
Both are smooth functions on T a . Of course the function G C°°(M\r) is continuous 
on M if and only if for all indices a 

Va; G T a : <j>a( x ) = 4>a( x )- 
Similarly, G C 1 (M) if and only if for all a 

\/x G T a : (j)+(x) = <j)~(x), and 
Vx G r ft : (n a 0) + (a;) = (n a (j))~(x). 

Here, (n^) 111 denotes the right and left limits of the normal derivative of 0. 

Let us define D a cf) to be the function on each manifold r Q given by 0+ — 0^ and let 
D na (p be the function given by (n a (f>) + — (n a (p)~ . Let Dcp and D n cf) be the corresponding 
functions in L°°(r). As the boundary of V has zero measure in V these functions are 
well defined and they are smooth on each T a . 

Let s,t < 0. The functional -F s r t (</>) defined by 

(i) FiM--=(\\mW) + \\DnH 2 Ht{r) ) l/2 , 

measures the continuity of the function and its normal derivatives in different Sobolev 
norms. Here the H s norm on V is defined as 

ll/ll^ s (r) = ll/all^ s (r a )- 

a 

Let now as above G C°°(M\f) and let x e C°°(M\T) be defined by = 
(A — \)4>(x) for x G M\r. Suppose that / G C°°(M) is any test function. Then, by 
Green's formula (see e.g. Prop. 4.1 in |Tay96|), 

/ </>{(A-\)f)(z)dii{z) = 

J M 

= [ ((A - X)(f))(x)f{x)dfi{x) + I (n(j){x)f(x) - <f>(x)nf(x)) dv{x) = 

JAf JdM\T 

= I X(x)f{x)dfi{x) + / (D n (/))(x)f(x)du r (x) - / (D^)(x)(n/)(x)^ r (i), 

where n denotes the unit normal vector field along T, dfi(x) = ^J~gdx is the measure 
induced by the Riemannian metric on M, and du(x), dur(x) are the measures on dM\T 
and T induced by the Riemannian metrics there. 
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If we think of as a distribution on M, then the above simply means that we have 
in the distributional sense 



(A - A)0 = X + D<f> S' T + D n (f> 5 T , 

where Sr is the Dirac delta distribution on T and S' r = nSr is its normal derivative so 
that for a test function / G C°°(M) the above distributions are given by 

{D n <j> S T ) (/) = J(D n <i>)(x)f(x)dvr{x) 

and 



(D<j>(S' T )) (f) = - jW)(x)(nf)(x)dMx). 
Note that the distribution 

h = D<j>6 r + D n (j) S r 

is in H~ 2 (M). Denote by F r (<f)) its if _2 -norm, where the iP-norm is defined by 

||^||^ ( m) = ||(1 + A) s/ ^|| 2 . 

By the trace theorem for Sobolev spaces (see e.g. Prop. 4.5 in |Tay96| ) the maps 
/ i — y f\r and / (n/)| r extend to bounded maps from H 2 (M) — > H 3 ^ 2 (T) and 
H 2 (M) ->■ H 1/2 (r) respectively. Hence, by duality, 

ms^H-^Kcmw^^ 

and 

||Ai0 8t\\h-*(M) < C'\\D n (f)\\ H _3 

for some constants C, C > 0. Thus, there exists a constant C > depending only on 
the geometry of M and T such that 

(2) ^(0)<^-^I_3(0). 

2 ' 2 

A constant C in this estimate can be explicitly found for a given manifold, for example 
by using gluing functions and the Fourier transform. We will illustrate in section [4] how 
to get a bound on C for a given hyperbolic surface when F consists of geodesic segments 
(cf. Q and (§). 

The method of particular solutions relies on the following simple observation. 



Theorem 2.1. Suppose that <p G C°°(M\r) C L 2 (M) with ||0|| 2 = 1 sitc/i */ia* 

Vx G M\T : (A - A)0(x) = xO), 
X GL 2 (M), < Hxll < »7- 
If F r ((j)) < e < 1, then there is an eigenvalue of A in the interval 

r^^' A+ i- e J - 
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Proof. As before let 

h = D(f)5' T + D n (j> 5 r 
Then g = (A + l)~ 1 h G L 2 (M) and by our assumptions 



Now note that 



Of course 



II ^ II 2 < c 

(A-A)(0- 5 )=x + (l + %. 
|x+(l + %|| 2 <^ + (l + A)|b|| 2 



and 



110 - 0II2 > 1 - \\g\U 

as ||0||2 = 1- From this it immediately follows that in case (A — A) -1 exists its operator 
norm is bounded from below by 

1 - IMI2 
77 +(1 + A) lulls' 

This implies the theorem since the operator norm of the resolvent at the point A is 
bounded from above by the inverse of the distance of A to the spectrum. □ 

Remark 2.2. The statement of the previous theorem also applies to situations where 
eigenvalues might be close to each other or have high multiplicities. Not to overload 
notation we state this here only as a remark. If there is an orthonormal set (0i)»=i... k 



in L (M) such that for each of the <f> the conditions of Theorem 2.1 hold, then the same 
proof shows that there are at least k eigenvalues (counting multiplicities) in a small 
interval around A. 

To apply the method of particular solutions one constructs normalized functions (f) 
which are eigenfunctions of the Laplace operator on each of the Mj with eigenvalue A 

and for which F r L 3(0) is small. The above theorems then show that A must be close 

2 • 2 

to an eigenvalue and <fi is close to an eigenf unction. We will show in the following how 
this can be done for oriented hyperbolic surfaces of genus g. 

3. Hyperbolic surfaces 

A hyperbolic surface is a 2-dimensional orientable Riemannian manifold endowed with 
a metric of constant negative curvature equal to —1. Such surfaces can be obtained by 
factorizing the hyperbolic plane EI = {x + iy G C | y > 0} with metric y~ 2 (dx 2 + dy 2 ) 
by a discrete co-compact hyperbolic subgroup T of SL(2, R). The isometric action of 
SX(2,R) is given by fractional linear transformations 

a b\ az + b 
c d cz + d' 







Another way to obtain a surface of genus g with hyperbolic metric is to glue 2g — 2 
hyperbolic pair of pants. A hyperbolic pair of pants is a genus surface with boundary 
S 1 OS 1 US 1 equipped with a metric of constant negative curvature —1 such that the 
boundary curves are geodesies (see Figure [T]). We will also refer to such a surface as a 
F-piece. 
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Figure 1. F-piece with boundary geodesies 

Such hyperbolic pairs of pants are, up to isometry, uniquely determined by the length 
of their boundary geodesies (^1,^2,^3)- Any hyperbolic surface can be decomposed into 
pairs of pants by cutting along 3g — 3 non-intersecting simple geodesies on the surface 
(see Figure [2]). This results in 2g — 2 pairs of pants. 




Figure 2. Pants decomposition of a surface of genus three into four F-pieces 

A pair of pants can be glued from a right angled geodesic octagon in the upper half 
space (see Figure |3| by identifying the sides b and h, as well as e and g respectively. 
As indicated in the figure the geodesic octagon is constructed from two right angled 
geodesic hexagons. 



B 




Figure 3. F-piece glued from a right angled geodesic octagon 
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Note that identification of the sides b and h only yields a subset of a hyperbolic 
cylinder so that every pair of pants can be constructed from a hyperbolic cylinder by 
cutting and gluing. 



4. The method of particular solutions on a surface of genus g 

4.1. Explicit error estimates for hyperbolic surfaces. The first ingredient in our 
algorithm is an estimate on the constant in the estimate ^ for F r in the particular case 
when M is a hyperbolic surface and T is a finite union of geodesic segments. To start 
let us assume that T consists of exactly one geodesic segment. The general case of a 
finite union of geodesic segments will later be reduced to this case. We can assume that 
M is realized as a quotient of the upper half plane EI by some co-compact hyperbolic 
subgroup in SL(2, R) and that D C EI is an (open) fundamental domain whose boundary 
is a geodesic polygon. This means we can identify functions on M with functions on EI 
that are invariant under the action of the group. We can also assume that the geodesic 
segment T is a segment on the imaginary axis. Then there is a tesselation of EI by 
translates Di of D. Let L > and suppose that D\, . . . , Djy is a finite number of 
translates of D such that dist(x, D) < L implies that x is contained in the closure of 
Vjf =1 Di. If L is chosen as the smallest distance between non-adjacent sides of D then 
N — 1 would be the minimal number of fundamental domains of D needed to cover an 
open neighborhood of D (and thus the number of neighboring fundamental domains 
sharing either an edge or a corner with D). Let \ : EI — > [0, 1] be smooth compactly 
supported function which is equal to one in an open neighborhood O of D with support 
contained in D = Vjf =l Di. Now suppose that h is a distribution on M. Then, 

\Kf)\ 

0^/ 6 C*°°(Af) || (A + l)/||i2(M) 



\Hh-*{M) = SUp 



sup IMX/)I 
0^/eC°°(Af) II (^ + l)x/IU 2 (D) 



Using 



(A + l)( x /) = /(A X ) + x(A + 1)/ + 2(V X , V/), 
where the inner product is with respect to the metric, we get 

|| (A + l) X f\\ LH D) < llx(A + l)/|| z2(6) + ||/A X || i2( o) + 2\\ (V X , V/) \\ L2(£)) . 
Of course 

2||(Vx,V/)|| i2(5) < 2||Vx|U||V/|| ia(fi) = 
= 2||Vx||oo V /(/,A/) L2(Z)) < ||Vx|UI|(A + l)f\\ LH6 y 



n 



Since / is invariant under the group action 



H(A + l)/llL 2( D) = v / iV||(A + l)/|| L2(D) , 

and therefore we obtain 

||(A + l) x /|| i2( ^ ) <C'||(A + l) x /|| L2(D) 

with 



C = V7V(l + ||Vx||oo + ||A X ||oo) 



Consequently, 



(3) \\h\\ H -HM) < C sup |/,(V/)| 



/eC°°(M) II(A + 1)(X/)II L2(5 ) 

To continue we use a particular coordinate system on the upper half space. Namely, 
we identify the upper half space with (— 7r/2,7r/2) x R using coordinates (tp,t) and the 
metric (cos tp)~ 2 (dtp 2 + dt 2 ) such that tp = coincides with the imaginary axis and 
the coordinate system is centered at the point i. Note that (tp, t) is related to Fermi 
coordinates (p,t) by (cosh p, t) = ((cos tp) -1 , t). We assume that in these coordinates T 
is identified with the segment {0} x [0,£] so that h has the form 

h(f) = h(f) + h 2 (f) = J F 1 (t)f(0,t)dt - J F 2 (t) (J^f) (0,t)dt. 

Moreover in these coordinates A has the form (cos 2 tp) A e where 

( d 2 d 2 
A e = - — + 



^ dp 2 dt 2 / 

is the Euclidean Laplace operator in coordinates (tp,t). Then we have 

||(A + l)(x/)fe (/)) = 

= |||c0S^|A e (x/)||l2 (( _ 7r/2i7r/2)xM) +2{xf, A e (x/))L2((- 7 r/2,7r/2)xffi) + 

+|||cos^r 1 x/||i2 (( _ 7r/2i7r/2)xK) > (C\\(A e + 1)(x/)||l2((-V2,V2)xr)) 2 , 
where C = inf{| cos tp\ \ xiv^) 0}- Collecting all terms we obtain 

H-3(M) < C G sup . 

/eC£°(]R 2 ) IK A e + 1J/||L2(R2) 
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However, 



Mf)\ 



1 



Mf)\ 



2vr j&i + f + v 

< 2ll^illi? 3 / 2 (R) 

\F 2 m\ 



2tJ r2 l + e + ^ 2 
1 



ll(Ae + l)/||, 

;i+e 2 +^ 2 )/(e^K^< 



ll(A. + l)/||, 



where and F 2 are understood as functions on K extended by zero from [0, £]. There- 
fore, we finally obtain 



(4) 



\h\\H- 2 (M) < C 



| Fl H2 



H 3 / 2 ( 



I II p 112 

ll-^llfl-VaflR)' 



where the constant C is given by 

1 



(5) 



C=-=( sup coshp)v / iV(l+||Vx||oo + ||A X || c 

(p,t)esuppx 



If T consists of several geodesic segments then of course the same inequality holds with 
the constant being the maximum of the constants for the individual segments. 

For low lying eigenvalues it is in practice easier to compute the L 2 -norms of F\ and 
F 2 . Of course since the negative Sobolev norms are dominated by the L 2 -norms the 
above gives estimates of ||/i||fr-2(Af) in terms of the L 2 -norms of F\ and F 2 . However, 



in this case a slightly better constant can be obtained using Theorems B.4 and B.5 
Indeed, equation [3] means by duality that 

||^||if- 2 (M) — C\\(A a + 1)~ /i||i2( H ), 



where Ae is the Laplace operator on the hyperbolic plane. By Theorem |B.4| and |B.5 
we have 



(6) 



||(A H + ly'hW^ < V /C 1 2 + C 2 V /||F 1 || 2 :2(M) + ||F 2 || 2 :2( 



where C% and C 2 are the constants in these theorems. 

Example 4.1. For a surface of genus 2 we use a right angled 12-gon and thus have 
N = 25. Let L be the minimum of the distance between non-adjacent sides. For each 
side Si we use a cutoff function Xi{p,t) = P(p/L), where P is the function 



P(x) 



1, x < 

0, x > 1 

l-3x 2 + 2x 3 , x G [0, 1], 
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and p are Fermi- coordinates with respect to the infinite extension of the side Sj oriented 
in such a way that the region p > does not intersect the polygon. Note that ||Vxi|| — ^ 
and Axi < j? + i^- With x — Xi ' ' ' X12 we then get 

(7) HVxIloo < §, 

(8) l|Ax||oc<g + |, 

where we used the well known formula A(xiX2) = A(xi)X2 + 2(Vxi, VX2) + A(% 2 )Xi 
and the fact that near a fixed point all but at most two of the functions x% are constant. 
Since x is not smooth we can not directly use it in our estimate. However, we can use 
a suitable regularization Xe and the fact the the second distributional derivatives of x 
are in L°° to show that the estimate above still holds with the function x- Collecting all 
terms we get 

WHh-hm) < 5V {CI + CI) + - + lj JWvw + WFtfwy 
Since 5^ (Cf + C|) < 12 this further simplifies into the easy to use estimate 

(9) WHh-hm) < 12 (j± + \ + 1) slW^^W^y 

This estimate improves by a factor of 4.69 if one is willing to believe the numerical 
values computed at the end of Appendix |Bj 

4.2. Basis functions on hyperbolic cylinders. Let I > 0. Then the hyperbolic 
cylinder Zi is the non-compact manifold obtained by factorizing the upper half space 

fe £ / 2 \ 

by the subgroup of SL(2, M) generated by the element . We will use Fermi 

y e~ e/2 J 

coordinates (p, t) which are related to the usual coordinates by 

(x,y) = e*(tanhp, — — ). 

cosh p 

For L > denote by the finite hyperbolic cylinder 

Z\ = {x E Z e I -L < p(x) < L}. 
For A E R let V W (Z^) be the space of functions / in C°°{Zf) that satisfy 

(A-A)/ = 0. 

Note that / G V {x) (Zf) has a Fourier expansion of the form 

f(P, t) = ( ak fa(P> *) + b k 1>k(ft, *)) , 

where {<pk)kez satisfy 

(A - A)0 fc = 
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and have initial data 



and the functions ipk satisfy 



and have initial data 



fc (t,O)=e 2 ^, 

d P Mt,o) = 0, 

(A - A)^ fc = 



Mt,o) = 0, 

These functions are of the form $fc(p)e 2?rlfc * //£ , where $fc(p) solves the ordinary differential 
equation 

. , , 1 d _ d Air 2 k 2 

10 — -coshp- + -J- - A)$ fc (p) = 

cosh p dp dp £ 2 cosh p 

with the corresponding initial conditions. A fundamental system of (non-normalized) 
solutions of this equation, consisting of an even and an odd function, can be given 
explicitly in terms of hypergeometric functions 

^„<=r,/ \ / , x 27rifc _ ,s nik 1 — s nik 1 , 9 . 
n ven (p) = (coshp)" 2Fl (- + —, — + —;-; - sinh 2 p), 

, , , . 2nik „ ,1 + s nik 2 — s nik 3 l9 . 

$f d (p) =smhp(coshp)— 2 F!(— + —, — + —;-; -smh 2 p), 

where A = s(l — s) (see |Borl0j . where these functions are analysed). Normalization 
gives the corresponding solutions to the initial value problems. 

Denote by be AN + 2 dimensional subspace of spanned by 4>k and ip k with 
\k\ < N. Given / e we can truncate the Fourier expansion to obtain an element 



f(N) e V W given by 



f (N) (p,t)= ^ {a k <f> k {p,t) + hMP,t))- 

\k\<N 



The C 1 -norm on the cylinder Z\ is defined as 

ll/ll^ = ll/lli- + ll^^/lli- + ll»p/llx- 

The following result is crucial for our approximation of eigenfunctions. 

Theorem 4.2. Suppose E V^ x \Zi) is bounded on Zg. Let L > and suppose that N 
is an integer such that ^osSl > ^' Then, 

U-<P [n) \\l^)<Un)U\\l 

u-^ N) \y {z -)<Px(N)\\<i>\\L~(z e) , 
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where X (N) = A 1 {N), (3 X (N) = a/ A^N) 2 + A 2 (N) 2 + A 3 (N) 2 , 

oo 

A 1 (N)=A (coshc^o))- 1 , 



jfc=iV+l 
oo 



2nk 



A 2 (N) = A ^(coshcfc^o))- 1 



k=N+l 



A 3 (N) = 4 



! Aix 2 k 2 



A 



fc=AT+l 



£ 2 cos 2 ((p ) 



(sinhc fc (^o)) 1 



and 



2nk 



2ttA: 



arccos 



sm(<yj) 



— vA arccot 



47T 2 fc 2 

P 



-A 



7^cos%)-~7 



Remark 4.3. 5oi/j f3\(N) and fl\(N) are decreasing exponentially fast in N for fixed 
A and groe very good bounds for reasonably large N (see Fig. [^p. 



Figure 4. /3\(N) on a logarithmic scale for A = 30, £ = 2 arccosh(3 + 
2v^2) and L = 3/2. The corresponding cylinder covers a fundamental 
domain for the Bolza surface. 



Proof of Theorem \4-£\ By rescaling we can assume in the proof that H^Hl 00 ^) < 1- Of 
course ip = <p — (j)( N > has a convergent expansion into a Fourier series 

(11) ip(p,t) = Y (a k (fi k (p,t) + b k ?p k (p,t)) 

\k\>N 

which converges uniformly on compact sets. Both <\> k and ip k are of the form Q k (p)e 2mkt / 1 
where solves the ordinary differential equation (10) with initial conditions either 
$ fc (0) = l,$' fe (0) = or $a,.(0) = O^'^O) = 1 respectively. Thus, $ fc is either even or 
odd with respect to the reflection p \-t —p. 
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Substituting <p = 2 arctan(tanhp/2) we write = ^^(v 9 )- Let 

4ir 2 k 2 A 



£ 2 cos 2 ip 



and define Sk = arccos |^ so that V& is non-negative on the interval [0, s k \- Then, by 
Lemma A. 2 we have for all k > N and <fi < s k : 

^k(sk) > ^k(v) coshc fc (v?) + %(ip)—. = sinh c k (<p) 

where 



^(v 9 ) = / y/V k {r)dr. 



This integral can be computed explicitly and gives the expression in the theorem. There- 
fore 

|*fc(¥>)| < (coshc fc (^))" 1 *fc(s Jfc ), 
\%(<p)\ < y/V k (ip) (sinhc fc (y>))~ 1 * fc (a fc ). 

Since the i/i(p,t) is bounded by 1 we have |afc||0fc(p, t)| < 1 and |6fc||^fc(p,i)| < 1. For 
any p > 0. With y?o = 2 arctan tanh ^ the above bound then implies that for any 
p G [— L, L] both |a/c||0fc(p,t)| and \bk\\ipk{p,t)\ are smaller equal than 

(cosh c fc (v?o)) _1 

and \a k \\^-(j)k{pi t)\ and \b k \\-^ip k (p,t)\ are smaller equal than 



\A4(<po) (sinhc fc (y? )) . 

Summing the Fourier series one obtains on Z L the bounds < Al, (JU/'I — ^2 and 
1^1 < As- □ 

4.3. Main estimates for the algorithm for a surface of genus g. Suppose that M 
is a surface of genus g. Then M can be decomposed into X and Y pieces (pairs of pants). 
Each of these pieces can be cut open along a shortest geodesic that connects two different 
boundary components to obtain a surface that can be identified with a closed subset of a 
hyperbolic cylinder (see Fig. [5J. In this way we obtain a collection (Mi)i = i r __ >m of closed 
subsets of hyperbolic cylinders Zg. with geodesic boundary components in such a way 
that the defining unique simple closed geodesic 7^ on is contained in Mj. Note that 
we use disjoint cylinders here, i.e. we use a different copy of Z^ even if the geodesies 
have the same length and correspond to the same model in the upper half space. As 
7i is a simple closed geodesic in M the subgroup it generates in the fundamental group 
will define a covering map Z^ — > M. Therefore, every eigenfunction on M lifts to a 
bounded smooth function on Ze i which is an eigenfunction of the Laplace operator on 
Z t: . 
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Figure 5. A F-piece that was cut open along two geodesic segments 
covered by a hyperbolic cylinder. The hyperbolic cylinder is obtained 
from the shaded region by identification of the two dashed circles. The 
inscribed finite cylinder is double shaded. 



The surface M can of course be constructed from the Mj by glueing, i.e. by identi- 
fying the different boundary geodesies. Suppose that we have a collection (T a ) a=li ^ M 
of geodesic segments and a collection (T & )a=i,...,M so that the geodesic segment T a is 
identified with the segment T & . Thus, T = can be thought of as a finite union 
of geodesic segments on M and M\T is the interior of the disjoint union JJ^ Mi of the 
Mj. This is exactly the setting described in section [2] where M\T gets identified with 
Y[ i Mi and the segments r a and T & form the boundary of M\T. Thus, for a function 
/ G C°° (LL Mi) which is smooth up to the boundary the boundary values at T a and 
Tq, represent the two different boundary values it has along T a if we think of it as a 
function on M. 

On each segment T a we choose a set {x ai i, x Qj2 , . . . , x a ^ na ) consisting of an even num- 
ber of equidistant points such that x aj i is the origin of the segment and x a ,n a is the 
endpoint. Let {1^1,1^2, . . . ,x &! n & } be the corresponding points on T & - We denote by 
S a the distance between neighbouring points on the segment r a and by 6 the maximum 
of (5 a ) a . 

Let N > be a fixed integer and suppose A > 0. Of course on each cylinder Zg,, and 
for each integer m such that \m\ < N there exists two linearly independent functions 
on this cylinder which are both eigenfunctions of the Laplace operator with eigenvalue 
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A and which are of the form 

2i7rm£^ 

= $(pi)e~5~, 

where are Fermi coordinates on the corresponding cylinder. We are choosing 

the basis functions orthonormal on the largest subset of Mj which is invariant under 
the rotation symmetry of the cylinder (for a F-piece this will be one half of a cut-off 
hyperbolic cylinder and corresponds to the double shaded region in Fig. [5]). Since such 
a subset will always have the form [L 1; L 2 ] x S 1 this can be achieved by choosing 
orthonormal on this interval. 

We extend these functions by zero to the disjoint union of all the cylinders and think 
of the collection of all these functions as a basis (^k)k=i,...,(4:N+2)m m the space 

We will use this basis to identify vectors in Q( 4N + 2 ) m with functions on the disjoint 
union of the cylinders, but also with L°° functions on M, by restricting the functions 
to the pieces Mj. For v E ^ m+2 ) m let ^(v) be the corresponding function. 
For each geodesic segment T a we form the matrices 



(Ai )Q! ) ik 
(A 2jQ )jfc ~ 
(A-3 t a)ik - 

where n is the outward normal unit vector field along the geodesic segment, t the 
normalized tangential vector field, and 6j are the coefficients in the composite Simpson 
rule Q i.e. 61 = l,b 2 = 4, 63 = 2, . . . , 6 n «-i = ^,K a -i = 1- Next arrange Ai )Q and 
A 2jQ into one matrix A a = A l a © A 2iQ , and add also the tangential derivatives A^, = 
Ai iQ © A2, a ©A3 jQ . All the matrices involved map into the same space and we understand 
the direct sum here as a direct sum of linear maps mapping into £.(' lN + 2 ) m _ Thus, the 
matrix Ai © A 2 is obtained from by attaching all the rows of A 2 as rows to the matrix 
Ai (similar to the Matlab notation [Ai; A 2 ]). 



i other coefficients such as those of the Gauss rule may also be used here and lead to analogous error 
estimates. 



5 a bi 



*fc(Xi), 



S a bi 



S n bi 
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Similarly, but with one change of sign 



*-2,a ik 



—^-n & W k (xi) 



{M,&)ik = Y -y-tfi*fe(a?i), 

and again A 5 = A 1; a © A 2 ,5, A~ = A^a © A 2 ,a © A 3) g. Now define 

A — $ tt A ai A = ©q,A q 
A — ®qAq, A — ©q,A^ 

as well as B A jv = A* — A*, B° N = A — A and Ca,at = A © A. Thus, Ba,at as compared 
to B° N contains also the tangential derivatives. 

For each v G (£(4JV+2)m vec ^ or go ^ v con tains the discretization of the jump 
of boundary data of ^(v) across T whereas Ca,atv contains the discretization of the 
boundary data itself. We will see below that the L 2 norm of an eigenfunction on M can 
be bounded from above and below by the L 2 -norm of the boundary data and thus the 
norm of Ca,atv will serve as a normalization in our method. To be more precise, the 
matrices are chosen such that the norms ||B° N v\\ and ||Ca,atv|| are the approximations 
of i^ (^(v)) and of the L 2 norm of the boundary data of ^(v) using the Simpson rule. 
The error of composite Simpson integration of a function g on an interval [a, b] with 
step size 5 is bounded from above by S |g~ a ^ supg 6 r a6 i |<7^(0I ( see e -§- |PT96j . p. 133). 
Thus, we have 

(12) ll|BVH 2 -Ko(*(v))) 2 |< 

< ^ (llA 4 (^(v) 2 )|r||^ ( r) + \\Dtn T (*(v) 2 )\ r \\ Lao{r) ) , 

(13) II|Ca,^vv|| 2 - (\\(*(v)\ T \\h {T} + ||n*(v))| r || 2 L2(r) ) | < 

< ^ (ll^(*(v) 2 )||L-cr) + ||AV(^(v) 2 )| r ||^(r)) , 

where D t denotes the tangential derivative. 

Moreover, as we chose the basis functions orthonormal with respect to the L 2 -inner 
product of a subset of M, we have 

(14) ||*(v)|U2 (M) > ||v||. 

Assume now that A is an eigenvalue of the Laplace operator and suppose that (f> 
is an L 2 -normalized eigenfunction. We may truncate the Fourier expansion of this 
eigenfunction on each of the cylinders to obtain a function < - 7V - ) which then corresponds 
to a fixed vector v such that (j)^ = ^(v). We choose N large enough such that Theorem 
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|4.2| applies and from the fact that is bounded by an explicitly computable constant 



see Appendix, Corollary E.3) the following estimates are explicit 



(15) \\(4> {N) - 4>)\r\\L- { r) < (3x(N)U\\oo, 

(16) ||(n0W-n0)| r || L< » ( r) < Px(N)U\U 

The error of the composite Simpson integration of a function g on an interval [a, b] with 
step size 5 is also bounded by 5\b — a\ sup^ 6 [ o6 ] \g'(0\ as one can easily see by replacing 
g in each sub-interval [ } of length 26 by the constant function g(x k +i)- 

Therefor^, we obtain 

H|CA,iW|| 2 - (||0|r||i 2( r) + I|n0| r ||! 2( r)) I < 2/3 A (iV) 2 ||0||L + 
+(2S£ r ) (||A(0 2 )|rlU-(r) + II An r (0 2 )|r|U~ ( r)) 



From the bound ( 23 ) of the Appendix and the fact that (f) is L 2 -normalized we conclude 
that 



|2 , 1 1 1 1 1 1 2 



|r|| L 2(r) I" ll n 0lr|| L 2 {r) J >ci(A), 

where c% is explicitly computable and depends only on A and the geometry of M and 
T. Using the bound L°° and the C l bound in the Appendix of eigenfunctions we can 
compute constants N c and 5 C such that for all N > N c and 5 < 5 C we have 

IIC vll> Cl(A) 

||^A,iVV|| > — - — • 

Note that this is a very crude estimate and we have used only the linear error estimate in 
Simpson's rule. This was done so that we can use the explicit estimates that we proved in 
the appendix. The critical constants S c and iV c guarantee that the numerical integration 
using the Simpson rule together with the approximation of (f) by the truncated function 
together yield a relative error of not more than 50%. For the method it will be enough 
to show that || Ca,a^ v || is bounded from above so that this estimate is sufficient. 

Let cri(B° N ) be the smallest singular value of B° N and ui(Ba,jv, Ca,aO the smallest 
generalized singular value (see e.g. |GL96j or also [Bct07j in the context of MPS) of 
the pair (B\ 7 n, C\,n)- Here the generalized singular values of Ba,at G Mat(m2 x mi) 
with respect to Ca,tv £ Mat(m 3 x mi) are defined as the singular values of B^.at as an 
operator from C mi to C m2 where the Euclidean norm || • || on C" 11 is replaced by the 
norm || • \\c XN defined by ||v||c AJV := ||Ca,jvv||. Note that 

ai(B A jv,CA,Ar) = inf BvaV " 



2 We use this estimate instead of ( 13 1 because a bound on the first derivative of eigenfunctions is 



easier to obtain explicitly and we do not need good accuracy for this bound 
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If N is large enough, so that the assumptions of Theorem 4.2 are satisfied, we have 

2^2 



ci(B 



A.TVj ^A,/V, 



< 



Cl (A) 



ft 



in case A is an eigenvalue. This follows immediately from the fact that an exact eigen- 
function satisfies the boundary conditions. 

Let as usual A = s(l — s) with Im(s) > and define a = Re(s) and r = Im(s). 
Assume that A' is not an eigenvalue and let s' , a' and r' defined as above. For each i we 
choose a hyperbolic cylinder Z L with L large enough such that Mj C Z L and such that 
Mi is relatively compact in the interior of Z L , so that the distance between dZ L and Mi 
greater than a positive number d > 0. For the truncated eigenfunction 0^ consider 



the boundary data 0^ © n<p ( z' 011 dZ L . Using the resolvent kernel of the Laplace 
operator ^ on the hyperbolic cylinder (see Appendix D ) we can reconstruct 0^ from 
its boundary data 



4>( N \x) = / kf(x,x')n x r(f)^(x') — (j)z(x')T^x'kf(x,x')dx'. 



Now define 



Then 



kf,(x, x')xi x :(j)^ \x') — 4>z \x')n x ikj{x, x')dx' . 



<f ) ( N \x)-<j ) ( f ) (x) = 

d d 
— k~ (x, x') rv0^ ( x> ) ~~ { x ') n x' ~Qzk~ {x 1 x')dx'ds. 



Using Lemma D.2 we obtain 



U0 (A ° 



4>J> \\c\M) <\S~ s'\C A 



■W.s..' I \\4>z'\\h{dZ) + \\ n( Pz \\L 2 (dZ) 

Of course the L 2 -norm of the boundary data is bounded from above by the L 2 -norm of 
the boundary data of the exact eigenfunction so that we obtain the bound 



110 



(TV) 



|C!(M) 



< Is 



5'|CMWvoi(r)||0|| 



Theorem 4.4. Wii/i the notation above we have for all N > N c and 5 < 5 C 



2^2 
ci(A) 



/3a(AO||0|| c 



j'lCM^VvoKnii^n 



i/A is an eigenvalue with normalized eigenfunction 0. Here 



and 11011c 1 are bounded 



by explicit constants (Corollary E.3 and E.5). 



Moreover, by Theorem 2.1 and the inequalities (12) and (14) we have: 



this is also often referred to as the Green's function of the Helmholz operator 
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Theorem 4.5. Suppose that v E C( 4N+2 ) m % s a unit vector and e = ||B° N v\\ and let 



T = C(t+ 6 -^ (||A 4 (^(v) 2 )||^ ( r) + |p t 4 n r (*(v) 2 )|r||^(r))) < 1, 



where C = C^C\ + Cf are the constants in Sections- 1\ Then, there exists an eigen- 
value in the interval [A — r , A + T ]. In particular the estimate holds if e = 
o"i(B° N ) and v is a corresponding normalized singular vector. 

Remark 4.6. The term e + ^ (||D t 4 (^(v) 2 ) || ioo(r) + ||D t 4 n r (*(v) 2 )| r |U- (r) ) can be 
minimized as another singular value problem, using the bounds on the derivatives of the 



basis functions that follow from the differential equation together with the Lemmata A. 4 



and \A.l , This results in a matrix, whose singular value is precisely the error estimate. 



Actual computations show that in practice the major contribution of the error comes 
from the first term. 



Remark 4.7. An analog of Theorem 4-4 holds for the n-th relative singular values By^ 
with respect to Cy^. In this case the term \s — s'\ gets replaced by the maximum of 
{\sk — s'\ | k = l...,n}, where Sk correspond to n pairwise orthogonal eigenvalues. 
We do not state this here as we want to keep the exposition simple. In the same way 
Theorem \4-5\ applies to situations with multiplicities and eigenvalues that are close to 



each other (see remark 2.2) 



4.4. Description of the algorithm for a surface of genus g. Since the fundamental 



solutions of the differential equation (10) satisfy the bounds of Lemmata A. 4 and A.l 
their higher derivatives can also be bounded explicitly simply by using the differential 
equation. We use coordinates r = sinh p and t on the hyperbolic cylinders, so that the 
differential equation is of the form 

(17) (-(1 + ^(1 + r 2 )^ + ^ - (1 + r 2 )A) S fc (r) = 0. 

Using subdivision into smaller intervals and the Taylor expansion around the boundary 
points of these intervals yields an approximation $[ M ^ (r) to the exact solution by piece- 
wise defined polynomial splines. Using the bounds on the higher derivatives the error in 
the C 1 -norm can be explicitly bounded from above. Since the measure on the hyperbolic 
cylinder is equal to drdt in these coordinates, normalization of the basis functions can 
be done within this error, simply by integrating these piecewise polynomials. In this 
way, the matrices B^, Ba,at and C\,n can be computed with explicit error bounds. 
The generalized singular value decomposition can then be done using any method from 
linear algebra. Once a numerical generalized singular value decomposition has been 
obtained the error of the singular values can be estimated from above as follows. Sup- 
pose that UDV is a numerically obtained generalized singular value decomposition of 
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B with respect to C. Thus, D is diagonal and U is unitary within a certain preci- 
sion. Moreover, V is unitary with respect to the inner product induced by C, again 
within a certain precision. Thus, ||U*U — 1|| < 8 X and ||V*C*CV - C*C|| < S 2 , and 
|| B — UDV|| < 63, where these three numbers reflect the errors. All three numbers 
can be estimated from above using interval arithmetics for example by bounding the 
operator norm by the Hilbert-Schmidt norm. Since the matrices Uo = U(U*\J)~ 1 / 2 
and V = V(V*C*CV)- 1 / 2 (C*C) 1 / 2 are unitary in the respective inner products the 
generalized singular values of UoDV are exactly the diagonal entries of D. By the 
characterization of generalized singular values of B as the norm distance of B to the 
rank k operators it follows that the error of the generalized singular values of B bounded 
by 

IKB-UoDVoXcrc)- 1 / 2 !!. 

Using functional calculus a longer calculation shows that this is bounded by 



01(C) V.2 01(C) 3 

where 01(C) is the smallest singular value of C. Here we assumed that 5 2 < 01(C). 
Note that the smallest singular value of C can also be estimated from below using 
exactly the same procedure for the ordinary singular value decomposition of C. Thus, 
interval arithmetics can be applied to obtain a lower bound for the generalized singular 
values of B with respect to C. 

An interval I is then tested for eigenvalues by computing a lower bound for the 



value of 0i(JZ x ,N)Cx,n) f° r a discrete set of points x in I. By theorem 4.4 for large 
enough N the distance between the points can be chosen such that eigenvalues can 
occur only near points where the lower bound is small enough. For a fixed e > 
this algorithm then will normally find a discrete set of points jji such that eigenvalue 
can occur only in e-neighborhoods of these points. The presence of eigenvalues and 
their multiplicities can then be tested using the singular values of B° jAr . An upper 
bound for small singular values can easily be obtained by first finding the singular value 
decomposition numerically and then using the numerically obtained singular vectors v 
to compute ||B° jyv||. Since the derivatives of the basis functions can be bounded using 



the differential equation (10) and the bounds of Lemmata A. 4 and A.l the error terms 
in Theorem 4J5 can directly be estimated from above. Using interval arithmetics one 
can therefore obtain rigorous interval inclusions for the eigenvalues. 

Remark 4.8. Once the presence of a single eigenvalue is established in a small interval 



theorem 4-4 can a ^ so be used to establish an error bound by narrowing the interval in 



which the eigenvalue can be located. 



24 



The analysis of rigorous error bounds can be simplified considerably by replacing 
the exact basis functions with the piecewise defined polynomial basis functions them- 
selves. These are of course orthogonal on any cylinder, thus normalization involves only 
integration of piecewise defined polynomials which can be done explicitly. Thus, the 
quasi-mode becomes a finite linear combination of products of piecewise defined poly- 
nomials in r and Fourier modes in t. The so constructed function does not satisfy the 
equation (A — A)$ = but instead (A — A)$ = x, where due to the nature of the differ- 
ential equation (1 + r 2 )x is again constructed out of piecewise defined polynomials and 
Fourier modes. Its L 2 -norm can very easily be estimated from above by integrating over 
the smallest cylinders that contain Mj. The fourth derivatives in the error of Simpson's 



rule in theorem 4.5 can also be estimated more directly in this case as a bound can 



easily obtained from the coefficients of the polynomial. 

5. Implementation and examples 
We have implemented our method in three different programs. 

5.1. Non-rigorous implementations. The first two programs are Fortran programs. 
One works solely for genus 2 surfaces and takes as coordinates mw-Fenchel-Nielsen 
coordinates. The second program takes as an input a labeled graph (encoding the way 
the surface is glued from F-pieces) and the Fenchel-Nielsen coordinates. Both programs 
then follow the algorithm described in this article. The basis functions are computed for 



each value of A by numerically solving the differential equation (10) using the DLSODA 



routine from ODEPACK (see (Hin831 [Pet83] ). Normalization of the basis functions 



is achieved by numerical integration of the solution of (10). In order to compute the 
singular values cr 1 (B° 7V ) we use LAPACK routines. The generalized singular values 
o"i(Ba,at, Ca,tv) are computed using QR-decomposition. A search algorithm then looks 
for the small minima of (^(B^at, Ca,at) as a function of A. A small enough step size 
guarantees that no eigenvalues are missed in the process. The other singular values are 
also necessary in the search: a small second singular value of B° N implies that either 
the eigenvalue has multiplicity greater than one, or another eigenvalue is very close. 
Our method can not distinguish between multiplicities and close eigenvalues if they can 
not be separated within the given precision. Both programs achieve accuracies close 
to machine precision for the low lying eigenvalues when the surface is far enough from 
the boundary of Teichmuller space (so that the DLSODA algorithm can compute the 
basis functions accurately enough). The code for the genus two program is available to 



the scientific community under GPLv3 and can be obtained from http: //www-staff 



lboro . ac . uk/~maas3/hyperbolic-surf aces/hypermodes . html 



Another program was written in Mathematica. It uses Taylor's method on intervals 



of size of order l/y/X in order to solve the differential equation (10) and compute the 
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matrices B^.at and Ca,at with arbitrarily high precision. The linear algebra such as QR- 
decomposition and search for the minimum is done within Mathematica. The achieved 
accuracy depends on the computing time and the machine precision used. We were able 
to compute the first 70 digits of the first eigenvalue for the Bolza surface (see Section 



5.3) using this program. 



5.2. Rigorous numerical implementations. To demonstrate the power of our method 
we modified the genus two Fortran and Mathematica codes to give rigorous error esti- 
mates for the eigenvalues. As we described in the previous section we use the Taylor 



series method to solve the differential equation (10) within a given accuracy and then 



follow the analysis described in section 4.4 The corresponding eigenvalue inclusions 



can be considered rigorous if one carefully includes rounding errors in the analysis and 
is willing to trust the implementation of basic mathematical functions in Fortran and 
Mathematica. 

5.3. The Bolza surface. In order to demonstrate our method we consider a special 
surface of genus 2. The so-called Bolza surface has Fenchel-Nielsen mw-coordinates 

(£ut l] e 2 ,t 2 ;£ 3 ,t 3 ) = 

= (2arccosh(3 + 2y/2), -; 2arccosh(l + y/2), 0; 2arccosh(l + v^O). 

This surface is known to maximize the length of the systole in genus 2 and it also 
maximizes the order of the symmetry group (see e.g. |Jen84} ISch93l ISch94] ) . The first 
eigenvalue on the Bolza surface was estimated by Jenni |Jen84] remarkably accurately 
to be in the interval [3.83,3.85]. Aurich and Steiner |AS89j (see also |ASS88j )give the 
value 3.8388 which was obtained using the finite element method. Figure |5.3| shows a 
plot of the singular value after QR-decomposition (N = 60, 5 = 0.001) of the matrix 
A © B and projection onto the first direct summand. Note that at the minima the 
function is very small and order 10~ 12 . 

The search algorithm in the Fortran program finds the first eigenvalue within 8 digits 
of precision. Using multiple precision in Mathematica we obtained the value. 

Ai « 3.8388872588421995185866224504354645970819150157 

where we believe that all digits are correct. The rigorous implementation that uses very 
simple bounds in its current form gives a rigorous error bound of 10~ 6 (which can be 
improved by simply investing more computing time). 

We believe that Ai is important as we conjecture that this is the maximal value 
of the first eigenvalue of the Laplace operator for constant negative curvature —1 in 
genus 2. This conjecture is supported by our calculations in Teichmuller space (the 
corresponding analysis will be discussed elsewhere) and is in line with the findings in 
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Figure 6. Smallest singular value after QR-decomposition for the Bolza 
surface as a function of A. 

( |JLN + 05j ). where it is shown ^ that the maximal value of the first eigenvalue of the 
Laplace operator with respect to any metric fixing the volume is attained for a singular 
metric on the Bolza surface ( [.TLN+05j ). 

6. Calculation of the Spectral Zeta function and the Spectral 

Determinant 

The spectral zeta function Ca(s) is defined as the meromorphic continuation of the 
function 

oo 

1=1 

where (Ai) ie N are the eigenvalues repeated according to their multiplicities. The above 
sum converges for 9R(s) > 1. As usual a meromorphic continuation to the whole complex 
plane is constructed using the Mellin transform 

1 f°° 

C A ( S ) = f ^y o f- 1 (tr(e~ A *) - 1) dt. 

and the fact that the heat kernel has an asymptotic expansion, by splitting the above 
integral into an integral over [0, 1] and one over [1, oo). From the heat expansion it fol- 
lows that zero is not a pole, which allows one to define the zeta-regularized determinant 
det^(A) of the Laplace operator by 

logdet c (A) = -Ck(0). 

4 this proof relies on a lower bound for the eigenvalue of a mixed Dirichlet-Neumann problem that 
was obtained numerically from a finite element method 
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In the examples we will also consider the value of the zeta function at the point —1/2. 
In physics this value is often referred to as the Casimir energy or vacuum energy. 
For our purposes we use the following splitting 

Ca(«) = p^y QV 1 (tr(e- A *) - l) dt + jT V 1 (tr(e" At ) - l) dt) , 

where e > 0. 

The Selberg trace formula applied to the heat trace (see e.g. Equ. (193) in |M04j ) is 

(18) trfe-") = V ° 1(M)e " ! r — dr + fY e —^ M 

4vrt J cosh 2 (7rr) f^ Z ^v / 4^2sinl 



n=l 7 v " m 2 

where the second sum in the second term is over the set of primitive closed geodesies 7. 
Splitting the integral and using the standard analytic continuation gives for 3?(s) > — N 

CaOO = j^yCZT (a) + Z?» + 7?"( a ) + 7?%)), 

where 



i=i 



iV 



s+k—l 



s + k-1 

k=0 



ra=l 7 



= EE/ t- 1 ^ '™*' " eft, 



where 



I^r) = ^V 2 U-^+i)* - g ^(r 2 + J dt, 
the 0^ are given by 



Vol(M) r°(-l) fc 7r(r 2 + l/4) fc 
ak -^Th k\ cosher) 

and r(x,?/) denotes the incomplete Gamma function 

/•oo 

r(x,y) = / t x - l e- l dt. 

Jy 

Differentiation results in the following formula for the spectral determinant. 



- logdet ? A = Ca(0) = L x + L 2 + L 3 

28 



where 



X 



J™ sech 2 (7rr) ^ ^ 2 (e{r 2 + \)) + + log(e ( r 2 + 1/4))) j d r, 

L t=\-\- e -*/4 ^ « df 

3 i^i^Jo 4^ 3 / 2 sinh(i^( T )) ' 

and E 2 (x) is the generalized exponential integral which equals x T(— l,x). All the 
integrals have analytic integrands and can be truncated with exponentially small error. 
They can therefore be evaluated to high accuracy using numerical integration. 

For fixed s and e > not too small the sums over the eigenvalues converge rather 
quickly and therefore T±(s) and L\ can be computed quite accurately from the first 
eigenvalues only. The error made by summing only over finitely many eigenvalues can 



be estimated directly by Lemma E.2 For not too large e > the sums over the length 
spectrum are also very fast convergent. As they are exponentially decaying as e tends 
to zero one can either neglect this contribution by choosing e small enough or one can 
compute them using a finite part of the length spectrum. The error introduced in this 
way is also exponentially decreasing as e — > and it can be explicitly estimated using 
the Selberg trace formula as follows. If the length of the shortest closed geodesic is L 
the second term in the Selberg trace formula Equ. (18) is for allt < T < \/L 2 + 1 — 1 
bounded by 

(19) F T (t) = ^tr(e" AT )e? + ^e - ^ 



and the heat trace for a fixed T can be estimated by Lemma |E.2| Therefore, we have 
for all T > e 



\T!' N (s)\ < f t^Fr^dt 
Jo 

\Ll\ < £ ±F T (t)dt, 



and the right hand side converges exponentially fast to zero as e — > 0. An analogous 
formula holds for the error if elements in the length spectrum up to length L are taken 
into account. 

Of course the short geodesies dominate the contribution from the length spectrum. 
Note however that as a consequence of the collar theorem there are at most 3g — 3 closed 
geodesies of length smaller that 2arcsinh(l) w 1.7627 on a surface of genus g > 1 (see 
e.g. |Bus92| ). 



29 



7. Spectral Determinants and Casimir Energy for special surfaces 

In this section we compute the values of the Casimir energy and the spectral deter- 
minant for the three isolated surfaces of genus two with large symmetry group as well 
as for other examples of genus two and genus three surfaces that have been treated in 
the literature before. All Fenchel- Nielsen coordinates are given in mu>-form, i.e. the 
3-valent graph describing this decomposition consists of two vertices that are connected 
by three edges. For all the surfaces we first used the Fortran double precision imple- 
mentation of our program to generate a list of eigenvalues. Even though our algorithm 
allows us to choose the stepsize in such a way that we are guaranteed not to miss 
an eigenvalue the current implementation does not achieve that automatically. In our 
computations we have chosen to check completeness of the list of eigenvalues by using 
the Selberg trace formula applied to the heat trace. This is possible as the absolute 
value of the contribution of the length spectrum in Selberg's trace formula is bounded 



from above by Fx(t), defined in Equ. (19), so that missed eigenvalues can easily be 



detected comparing the heat trace computed from the Selberg trace formula. Note that 
unlike heuristic checks using Weyl's law this gives a rigorous way to check complete- 
ness up to a certain value. The data files for all the examples below can be found at 



http : //www- staff . lboro . ac .uk/~maas3/publications/eigdata/dataf ile .html 



7.1. The Bolza surface. The Bolza surface has symmetry group 5*4 x Z2. Its order, 
48, thus maximizes the order of the symmetry group amongst all genus 2 hyperbolic 
surfaces. The Bolza surface is referred to as the regular octagon (and sometimes as 
the Hadamard-Gutzwiller model) in the physics literature as it can be obtained from a 
regular octagon in hyperbolic space by identifying opposite sides. Its Fenchel-Nielsen 
coordinates are 

(*i,ti;4i,*2;4»,*3) = 

= (2arccosh(3 + 2y/2), -; 2arccosh(l + V2), 0; 2arccosh(l + v^),0). 

We computed the spectral determinant as well as the Casimir energy ((—1/2) and 
obtained 

det c (A) « 4.72273280444557, 
Ca(-1/2) « -0.65000636917383, 

where we believe that all decimal places are correct. This accuracy was obtained by 
computing the first 10 eigenvalues with a 30 digit precision, the first 500 eigenvalues 
in quad precision and subsequently by using the simple length spectrum up to £ = 9. 
We relied on the multiplicities of the length spectrum as computed in |AS88j . Note 
that the spectral determinant and the Casimir energy can both be computed accurately 
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within 5 decimal places without taking the length spectrum into account at all, using 
the eigenvalues only. Completeness of the list of the first 500 eigenvalues was checked 
using a list of 700 computed eigenvalues. To compute the first 500 eigenvalues of the 
Bolza surface to a precision of 12 digits about 10000 A-evaluations of generalized singular 
value decomposition were needed. This took about 10 minutes on a 2.5 GHz Intel Core 
i5 quad core processor (where parallelization was used). 



Figure 7. Ca(s) as a function of s for the Bolza surface 

We conjecture that the spectral determinant is maximized in genus 2 for the Bolza 
surface and we therefore believe that the above number is important. Our computations 
show that Bolza surface as a point in Teichmiiller space is indeed a local maximum. Note 
that for the Bolza surface is known to be a critical point for the spectral determinant. 

7.2. The surface with symmetry group D§ x This surface has Fenchel-Nielsen 
coordinates (see [ BSQ 5], sec. 3.5) 

{i\ ,h ; £ 2 M 4, h) = (2 arccosh(2) , 0; 2 arccosh(2) , 0; 2 arccosh(2) , 0) . 

It is also known to be a critical point for any modular invariant such as the Casimir 
energy or the spectral determinant. Computing the first 500 eigenvalues with double 
precision and neglecting the contribution from the length spectrum other than from the 
12 primitive closed geodesies of length 2 arccosh(2) we obtain 

det c (A) » 4.428000668, 
Ca(-1/2) « -0.67250924. 
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7.3. The surface with symmetry group Z 5 x Z 2 . This surface has Fenchel-Nielsen 

coordinates (see [.BSD.5], sec. 3.5) 

= (2 arccosh^, 0; 2 arccosh^, ±; 2 arccosh^±M, ± 
v 2 2 2 4 2 

It is also known to be a critical point for any modular invariant such as the Casimir 

energy or the spectral determinant. Computing the first 500 eigenvalues with double 

precision and neglecting the contribution from the length spectrum we obtain 

det c (A) « 4.358289, 
Ca(-1/2) « -0.67914508. 

7.4. The Gutzwiller octagon. This is a surface obtained from a regular hyperbolic 
octagon in the upper half space using a side identification different from that for the 
Bolza surface (see |Nin95] where also the low lying spectrum of this surface was com- 
puted). The Fenchel-Nielsen coordinates of this surface are given by 

(4,ti;4,*2;4,* 3 ) = 

= (4 arccosh^, I 2 arccosh^, 1; 2 arccosh^, h 
1 v/2 4 ^2 2' ^2 '2 ; 

Computing the first 500 eigenvalues with double precision and neglecting the contri- 
bution from the length spectrum we obtain 

det c (A) « 3.76048, 
Ca(-1/2) ~ -0.72747. 

7.5. The example of Aurich and Steiner. Eigenvalue statistics of a certain generic 
octagon was investigated in |AS93j . The Fenchel-Nielsen parameters of this surface are 
roughly (see |ABCKS| ) Q 

(i^h) « (3.717414183638,0.0304758025243), 
(£ 2 ,t 2 ) « (3.303402988815,0.2711895405183), 
(hM) ~ (3.408324727953,0.2187479424048). 

Computing the first 500 eigenvalues with double precision and neglecting the contribu- 
tion from the length spectrum we obtain 

det f (A) » 3.959168, 
Ca(-1/2) « -0.715195. 



"■see Acknowledgements 
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7.6. Examples of Rocha and Pollicott. These are the examples treated in |PR97j . 
where approximate values for their spectral determinants are given. 

The first example has Fenchel-Nielsen coordinates 

{£i,t 1 ;e 2 ,t 2 ;£ 3 ,t 3 ) = (5.05,0; 1,0; 0.9,0). 

Computing the first 500 eigenvalues with double precision and taking into acount only 
the short primitive geodesies in the length spectrum we obtain 

det f (A) « 0.395833, 
Ca(-1/2) « -1.817507. 

The second example has Fenchel-Nielsen coordinates 

{£ u h; £ 2 ,t 2 ; £ 3 , t 3 ) = (0.98, 0; 3.5, 0; 0.98, 0). 

Computing the first 500 eigenvalues with double precision and taking into acount only 
the short primitive geodesies in the length spectrum we obtain 

det c (A) « 0.6114618, 
Ca(-1/2) « -1.6541313. 

The third example has Fenchel-Nielsen coordinates 

{£xM\£2:t 2 ;£ 3 ,t 3 ) = (5,0; 1,0; 1,0). 

Computing the first 500 eigenvalues with double precision and taking into acount only 
the short primitive geodesies in the length spectrum we obtain 

det f (A) « 0.5124672, 
Ca(-1/2) « -1.6591527. 

The values do not agree with the values obtained in |PR97j . We believe that not 
enough lengths of primitive geodesies were taken into account there. 

7.7. An example of genus three. The following example has genus three and is 
a double cover of the Bolza surface. It is obtained by cutting the two copies of the 
Bolza surface open along the two geodesies of length 2 arccosh(l + a/2) used in its pants 
decomposition (see section 7.2 ) The two resulting surfaces of type (0, 4) are then glued 
by identifying each boundary with the the corresponding boundary on the other copy 
of the surface. We computed the first 600 eigenvalues on this surface and obtained 

det c (A) « 9.6507, 

Ca(-1/2) « -0.7802. 

For this accuracy the contribution of the length spectrum can be neglected. 
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Appendix A. Estimates on the basis functions 

This section deals with growth estimates and bounds on the function that sat- 
isfies the differential equation 

^ ^ ^ cosh pdp C ° S ^ 'dp cosh 2 p ^ ^ 

with initial conditions 

$(0) = a, -f$(0) = b. 
dp 

Substituting <p = 2 arctan(tanhp/2) this can also be written as 

/ d? , A 

Gt(/2 2 COS (/? 

where = ^(piv)) satisfies ^(0) = a and ^-^(0) = b. The derivatives in different 

coordinates are related by coshp^ = Now define a potential 

A 



V(<p) = A- 



cos 2 ip ' 



Lemma A.l. Suppose that A > and V(ip) does not vanish on an interval I C [0, 7r/2). 
Tnen £ne function E(<p) = \^/(ip)\ 2 — yj^\^/' (ip)\ 2 is non-increasing on I. 

Proof. Simple differentiation gives 

i , . 1 2Asin<z>, . , 9 

E V = _ 7t7TTi2 3"^ * V ^ °> 

COS J (/9 

which implies the statement. □ 

Lemma A. 2. Suppose that A > and V(ip) > on some interval [ipo,<pi] C [0, 7r/2) 
and assume $(</?o) > and $'(Vo) > 0- Define 

J fa 

Then for all ip G [</?o, we have 

$(<p) > Hfo) cosh (c Vo (<p)) + -j2== &((p ) sinh (cvo(^)) . 

and 



> a/^M^o) sinh ( Cw) ( V 9)) + y^^$ / (^o) cosh (c V0 M). 

Proof. Since every solution can be obtained as a sum of one with $(y?o) = and one 
with <J>'(v?o) — we can assume without loss of generality that either $(y? ) = or 
$'(</?o) = 0. We also assume that $(</?) is non-zero. Define 

1 



h((p) := $(y? ) cosh (c^oH) + / —— &(<po) sinh (c^H) • 
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Then h(ip) and $(y?) have the same initial data at ip . It follows from the differential 
equation and the assumption $(</?o) > and $'(^0) > that $(<^), $'(y?) an d ^"(v 9 ) are 
positive on some interval (</? , 5). Moreover, where <3>(y?) > and <&'((f) > the function 
$((/?) is convex and therefore, $(99), $'(<^) and $"(y?) are non-decreasing functions on 
[<Po, ifii]. In particular $(y?) > on (<£> , The differential equation implies 

The right hand side of this is non-negative as on can easily see by direct calculation 
that 

where T{ip) is either tanh((/?) or coth((^) depending on the inital conditions. Thus, 

(in = m<p)H<p) - h'teMv)) > 

on (</?o, <fii] and consequently 3>(v?) > h((p). The equation also implies that $(y?) — ft.(<^) 
is non-decreasing which gives the inequality for the derivatives. □ 

The following lemma is a direct consequence of Gronwall's inequality applied to the 
first order system 

A ( ) = ( c )( $(v?) 

Lemma A. 3. Suppose that A > 0, c > and \V(ip)\ < c 2 on some interval [ipo,(pi] C 
[0, 7r/2) Then for all ip G [<po, (pi\ we have 

(MV)\ 2 + < e c| ^ ol (|«f (V^o)l 2 + ^|$'M| 2 )i 

For intervals (<p , ipi) where V(ip) is positive this can be slightly refined. 

Lemma A. 4. Suppose that A > and V(ip) > on some interval (</?o,</?i) C [0, 7r/2) 
Then for all ip e (<po, <fi) we have 

Proof. The equation is equivalent to the system 

\ _( -/V(^j\ ( *(</>) 

The operator norm of the family of matrices in this system at the point ip is given by 




2\/(v?) V v 2y(^)- 
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which is bounded from below by 



given that V'(ip) < 0. The lemma now follows immediately from Gronwall's inequality. 

□ 

Appendix B. Estimates for the resolvent on the upper half space 

Let Ah be the Laplace operator on the upper half space equipped with the Poincare 
metric. Then the integral kernel k s (z,w) of the resolvent 

Rs(is)(A m ) = (A H -s(l -s))- 1 

is given by 

k s (z,w) = <5- s (cosh dist(z, w)). 



cosh dist(z,w) — 1 
2 



Here Q s is the Legendre Q-function. Using the point pair invariant u(z, w) 

k s (z, w) = ^-Q- S (l + 2u(z, w)) = FJu(z, w)), 
Ztt 

where F s (u) is given by an integral (sec e.g. [Iwa02j) 

F s (u) := ^g_ s (coshp) = 1- [ (£(1 - Or ^ + u)- s d£. 

This formula can be used to construct a meromorphic continuation of the resolvent as a 
function of s to the whole complex plane. We will need the following estimates on this 
function all of which follow right from the above integral formula. 

Lemma B.l. Suppose that a = Re(s) > 0. Then 

, d \s\ (r(q)) 2 _ a _ x 

— r 3 \U)\ < ; 

l du v ;| - 47T r(2cr) 

Lemma B.2. Suppose that a = Re(s) > 1. Then 

. . , . 1 1 + u 

\F 8 {u)\ < — log , 

47T U 

i d „ , s , \s\ 1 
-FJu) < 



du An u(l + u) 

A similar estimate as in the above lemma also holds for Re(s) > (c.f. |Iwa02j . 
Lemma 1.7, where however (1.48) seems to be inaccurate) 
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Lemma B.3. Suppose that a = Re(s) > 0. Then, there are constants Ci(cr), C^{a), 
depending only on a, and a constant C^s), depending only on s, such that 

|F»|<-^log— + C 1 (a), 
Aix u 

\^F s {u)\<^ l - + C^). 
du Air u 

Moreover, C^s) can be bounded by C^{a)\s\, where C^(a) depends only on a. 

In the important case a = 1/2 one can choose C*i(l/2) = 1, C*2(l/2) = 2 and 
63(1/2) = 1, although these choices are not optimal. The statement that C^s) can be 
bounded by C 4 (cr) |s| is not optimal either for large values of 3(s), where cancellations 
determine the asymptotic behaviour as u — > 0. 

Theorem B.4. Let T be a geodesic segment in the upper half space and let (p,t) be 
Fermi coordinates with respect to the infinite extension ofT. Suppose a distribution hi 
is defined as 

h(f)= [ F 1 (t)f(0,t)dt, 



where F x G L 2 (R). Then 

IKAh + s^)- 1 /^!!^^) < c 1 ||F 1 || L2( 



where the constant C\ is given by (21) and satisfies C\ < 1.75. 



Proof. Since the Laplace operator commutes with the SL(2, M)-action we can assume 
without loss of generality that the geodesic segment is part of the imaginary axis. The 
function (Ae + 3/2)~ 1 hi can be computed by convolving hi with the integral kernel 
k s (z, w) for s = 3/2. We use Fermi coordinates (p, t) so that the integral kernel is given 
by k s (p,t;p',t') = k s (p,t- t';p',0). Thus, 

((A H + 3/2)- 1 /i 1 )(p,t) = J F 1 (t')k s (p,t-t',0,0)dt'. 

Since this is a convolution the generalized Young inequality applies so that 

|((A e + 3/2)- 1 /i 1 )(p,t)| 2 rft< HFxIl^ (J \k e {p,t,0,0)\dt" 
The value of the point pair invariant u(z, w) at w = i in these coordinates is given by 

u(p, t, 0,0) = - (cosh t cosh p — 1) , 
so that the inequality follows with 

i 

(21) Ci= ( l{ I \F S (- (coshtcoshp- 1)) \dt) 2 cosh pdp 
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The integral over t can be estimated for cosh p > 3 as 

\F s (l (coshtcoshp - l))\dt < ^ (r( 5 Q 2 , ))2 (2) 3 / 2 (coshp)- 3 / 2 /°°(el*l/2 - l/3)~Ut = 
2 4tt 1 (3) J 

= - — 7r — 2 arcsin \ + 2 arctan (cosh p)~ 3 ^ 2 

For cosh p < 3 we split the integral over t into two parts. One with cosh t > 3 and one 
with cosht < 3. The former gives 



L 



(cosh t cosh p — < 



cosh t>3 

= ^ U/2 + log tanh Q log(3 + 2^)) ) 



and the latter 



-, \ \ i i 1/" , cosh t + 1 , 
|F S (- (coshtcoshp - l))\dt < — / log — -dt < 

cosht<3 2 47T J cos ht<3 COSh t — 1 

^ ;-arccosh3 ^2 , 4 arccosh(3) log ( 1 H T73)2 ) + 4 arctan (|arccosh(3)) 

< n / log ,., di V au.o,,,, , 



2tt J t 2 2tt 

Squaring and integrating finally gives 



Ci < (3 - 2v^) a? + v 7 ^) 2 « 1.7485475 

where 

ai = 3tt + 6 arccot ^2v^) - 12 V2, 
02 = ^ + log ^tanh Q log (3 + 2V^)^ + 
8 (arccosh(3) log (l + areco ^ h(3)2 j + 4 arctan (|arccosh(3))) 



7T 



□ 



Theorem B.5. Let V be a geodesic segment in the upper half space suppose (p,t) are 
Fermi coordinates with respect to the infinite extension ofT. Suppose a distribution hi 
is defined as 

h 2 (f) = - J o F 2 (t)(—f)(0,t)dt. 

Then 

|| (Ah + 1) 1 ^2||z / 2( H) < C 2 1 1 i*2 1 1 z, 2 (o,£) , 



where the constant C2 is given by (22) and satisfies C2 < 1.61 
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Proof. As in the proof of the previous theorem we can assume without loss of gener- 
ality that the geodesic segment is part of the imaginary axis. As before we use Fermi 
coordinates (p, t) so that the integral kernel is given by k s (p, t; p', t') = k s (p, t — t'; p', 0). 
Thus, 

(Ah + S^)- 1 ^)^^) = J F 1 {t'){-d p/ k s {p,t-t',p',0))\ p/=o dt'. 
Now by the chain rule 

{-d,k.(p, t - t', p', 0)) \ p/=0 = (-^F' s (u))(p, t - t', 0, 0) = s^F'Mp, t - t', 0, 0)). 

Again, by the generalized Young inequality 

|((AH + 3/2)-%)(p,t)| 2 rft=||F 1 ||i 2 (^) 2 (| \F' s (u(p,t,0,0))\dty . 



Thus, 

|| (Ah + 1) _1 ^2|||2( H ) < C 2 \\F 2 \\ L 2( ^, 



where 

(22) C 2 =(f(^J |^Kp,t,0,0))|^y(^) 2 cosh pc /p 

Splitting the integral into a region cosh p > 3 and cosh p < 3 and using the estimates 
IB.ll and lB.2l one obtains 

- 16 



^ (27 - 16v^) (sV2 - 6tt + 12 arctan ( v^) ) 2 + 48 (2V2 + arccosh(3) N 
The numerical value of the right hand side of this inequality is about 1.6086. □ 



Both constants C\ and C 2 can be computed using numerical integration and their 
values are about 

d « 0.313, C 2 « 0.343, 

where we believe that at least two digits are correct. The corresponding estimates are 
non-rigorous however. 

Appendix C. Estimates on derivatives of the resolvent 
We will need the following bounds on the derivative of k s with respect to s. 

Lemma C.l. Let F s (u) be the function defined in the previous section. Then, for any 
a > we have the estimates 

l^sHI < C^a) log(l + u)u-° + C 2 (*)u- ff , 



, ^q-/s{u)\ < C^u-*- 1 + |*| (Ci(a) log(l + u) + C 2 {aj) u°-\ 
\q^q-/s(u)\ < |2a + l|<7i(<7)u-'- 2 + \s(s + 1)| (C^a) log(l + u) + C 2 {aj) u^ 2 
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where 



and s = cr + ir. In particular (7i(l/2) = | and C 2 (l/2) = 
Proof. Simple differentiation shows that with s = a + ir: 

i^")i = is/iog^«(i-o)-« +ll) -^i< 

which yields the first estimate. The second estimate and the third estimate are obtained 
in a similar manner from 

d 2 



F s (n)\<^J\ai-or-\i+ur- i di 



and 



duds 

-g/ 1 i° g (f^)tt(i-or'K + »)-'-^. 
i^.(«)i<^/tt(i-or'tt + «)— « 

,s(s + 1)1 r io g (^f^)(«i - £))-'(? + 1 .)-"- 2 <i5. 
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Near zero one obtains similarly. 



□ 



Lemma C.2. Lei F s (-u) fre the function defined in the previous section. Then, for any 
a > we have the estimates 

\^- s F s (u)\<C 3 (a), 
\uJ^- s F s (u)\ < C 4 (a) + C 3 (o-)(r 2 + a 2 ) 1 / 2 , 

where s = a + ir. 

In the important case a = 1/2 one obtains 

^3(1/2) = J, 

C 4 (l/2) = sup — / - U — _ d£. 
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The latter constant can numerically be calculated and is approximately given by 

64(1/2) « 0.167878. 



Appendix D. Resolvent estimates on hyperbolic cylinders 



the subgroup of SL(2,M) generated by the element [ ^ _ ). We use Fermi- 



Let Z be the hyperbolic cylinder obtained by factorizing the upper half space by 

V/ 2 o 



coordinates (p,t) which identifies the cylinder with M. x R/(£Z). The integral kernel 
of the resolvent of the Laplace operator Az can be constructed from the resolvent of 
Ah by averaging over the group: Let x = (p, t) and x' = (p f , if) be distinct points on 
the hyperbolic cylinder and denote by k s (x, x') the integral kernel of (Ah — s(l — 
Then the integral kernel kf(x,x') of (Az — s(l — is given by 

kf((p, t),(p',t')) = Yl k *((P> *). V, t' + n£)) = Yl k ^P' 1 ~ ^ n£ )) 

The sum converges absolutely whenever dt(s) > because of the estimate in Lemma 



B.l and the fact that the distance between the points (p,t) and (p',f + n£) grows like 



nt as \n\ — > oo. 



Combining the estimates in Lemma C.I with Lemma B.3 we obtain 



Lemma D.l. For any compact subset M C Z we have on MxM the following estimates 

\kf(x,x')\ < ^\logp(x,x')\+Cf(a), 

r M (?) i 

| g rad^f(x,xO|<^li^ +C 3"W, 
where p(x,x') = 2 arcsinh(w(x, x') 1 ^ 2 ) is the distance between x and x' and a = 3?(s). 



Similarly Lemma C.l and Lemma C.2 imply the following estimates 



Lemma D.2. Let Z L C Z be a truncated hyperbolic cylinder and let d > 0. Then for 
all points x G Z L and x' G dZ L with p(x, x') > d the following estimates hold. 

\^kf(x,x')\<C^ d (a), 
|^grad x ,fcf (x,x')\ < C 2 M,d (s), 

d 

|— gTad x ,gmd x kf(x,x')\ < C^ I4 (s). 

Suppose now that M C Z is a subset of the hyperbolic cylinder Z whose boundary 
DM is a finite union of geodesic segments. If ip is a function on M, smooth up the the 
boundary such that 

(A z — s(l — s))if? = 
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in the interior of M. Let 0i and 02 be the restriction of ip to the boundary and its 
outward normal derivative respectively. Of course, as functions on Z we have 

in the sense of distributions, where h if the distribution defined as 



Kf)-= \ (<9 n /)(z)0i(z) + f{x)(j) 2 (x)dx. 

JdM 

From this it follows that ip m the interior of M can be expressed as 

i> = {& z -s{l-s))- 1 h. 
This is sometimes written in integral form 

ip(x) = / d n ykf (x, x')4>i{x') + kf (x, x')4>2{x')dx' . 

JdM 



Using, Lemma D.l and applying the generalized Young inequality we obtain the follow- 
ing estimate. 

(23) IHU 2 (M) < C M {s) ( \\(pi\\ 2 L 2 {dM ) 

L 2 (OM)j 

Appendix E. Appendix-Estimates on the Counting Function 
Let H be the Heaviside step function defined by 

{1 y > 
1/2 y = 
y<0 

Suppose that X is a compact connected oriented hyperbolic surface and let (0j) ig N o be 
an ortho normal basis of eigenfunctions for the Laplace operator on X. Let 

= A < Ai < A 2 < . . . 

be the corresponding sequence of eigenvalues. Let e(x,y,r) be integral kernel of the 
operator sign(r)H(r 2 — A) and denote by N x (r) = e(x,x,r) be its restriction to the 
diagonal. Of course, by the spectral theorem, 

N x (t) = sign^^H^-A,)!^^)! 2 - 

Integrating this over X yields the eigenvalue counting function N(t): 

iV(r) = sign(r)^H(r 2 -A l ). 

i 

By construction N x and iV are non-decreasing odd functions on R. Let d(x) be twice the 
injectivity radius at the point x E X. If L is the length of a systole we have the estimate 
r(x) < L. Using Selberg's pretrace trace formula and finite propagation speed one can 
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easily see that the cosine transform of N' x (r) coincides on the interval (— d(x), d(x)) with 
the cosine transform of the function F'(r), where F(t) is defined by 



F(0) = 0, 

F'(t) = 2^-E(t 2 - l/4)|t| tanli^^/i 2 - 1/4). 
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Note that since 

poo y 

/ 2t(tanh(7rv/t 2 - 1/4) - l)dt = 

and tanh(7TA/t 2 - 1/4) - 1 is non-positive on [1/2, oo) we obtain the estimate 

-~ < sign(t)G(t) < 0, 

where G(t) = F(t) — sign(t) ^t 2 . Let v be the first Dirichlet eigenvalue of the operator 
^ on the interval [—1/2, 1/2] and let (f>(x) be the corresponding normalized eigenfunc- 
tion. We think of as a function on K by extending by zero. It can be easily worked out 
that v is the smallest non-zero solution to the equation cosh(A) cos(A) = 1 for A > 0. 
Its numerical value is v w 4.73004074. The estimate v < 5 is an immediate consequence 
if the intermediate value theorem. As a test function in the Fourier Tauberian theorem 
we use the function p = (0) 2 . Let p$ and ps t o be defined as in [SafOl] as 

ps(r) = 5p(8r), 
Ps,o( r ) = <*Pi,o(<fr), 

where 

POO 

Pi,o(r)= / tp(t)dt. 

J r 

Using psflir) > and ||/0($||n = 1 we obtain 

(24) ( Ps *G)(r)< \\G\\oo = ^ 

(25) (ps,o*G')(r)<0. 

The Fourier Tauberian Theorem 1.3 in [SafOlj together with the estimates Equations 



(2.9) and (2.10) in [SafOl] and our estimates (24) and (25) imply the following bounds 
of N. 

Theorem E.l. The counting function on a hyperbolic surface satisfies. 
NJr) < — r 2 + 7r^{r + 



An \ nd(x) d(x) 3 

NJr) > — [r 2 — (r + 



4/T \ ird(x) d(x) 3 
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Define the function N X (X) = H(A)A^(y |A|). The local heat kernel trace k t (x) for 
t > may be defined as follows 



^\Mx)\ 2 . 



If this sum is cut off at a point c > the remainder is given by 



R$(x) = e" 



Aj U.fx)| 2 . 



Then of course integration by parts gives the following formula 



R c t (x) = - lim N x (c + e)e~ ct + t \ N x (X)e~ dX. 



Combining this with the estimate from theorem E.l we obtain 



Lemma E.2. The remainder R^(x) satisfies the estimate 

4nR c t (x) < e~ ct (b u {x) + B t {x) + (A u (x) + A l {x))^/c~ + ^ + 

+A u (x)^ t (l-eri(Vc~t)), 



where 



A u (x) 



4z/ 2 + 2vtt 



Mx) 



ird(x) 
4z/ 2 



ixd{x) 



B u (x) 

Bl{x) : 



4i/ 3 + 2v 2 ix 1 
+ 



Trd(x) 2 ' 3 
4z/ 3 1 
7rd(x) 2 3 



Setting c = one obtains a bound for the local heat kernel. The heat trace 



trie 



-At\ 



is obtained by integrating the local heat trace. Lemma E.2 therefore yields bounds on 
the heat trace in case c = and on 



R c t (x)dx. 



x 



in general. 



Theorem E.l also immediately gives bounds on the eigenfunctions 



Corollary E.3. Let ip(x) be an eigenfunction of the Laplace operator with eigenvalue 
X such that II^IIl 2 = 1- Then, 



i < 



(A 1/2 + — ) + 



47r \ ird(x) 



v 
d(x) 
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For numerical purposes we will use the bound 



where L is twice the radius of injectivity (which equals to the length of the shortest 
closed geodesic). 

Similarly, one gets a bound on the derivative of the eigenfunctions in the following 
way. Fix a point i6l and a unit vector n. As above we can define 

^ 1 (r)=sign(r)^H(r 2 -A J )|n0 l (x)| 2 . 

i 

Using the pre-trace formula and finite propagation speed it can easily worked out that 
cosine transform of (A^)'(r) coincides on the interval (— d(x), d(x)) with the cosine 
transform of the function F'(r), where F(t) is defined by 

F(0) = 0, 



F'(t) = — H(t 2 - l/4)|t| 3 tanh(7rv/t 2 - 1/4). 
Air 

Using 

17 



POO 

/ t 3 (tanh(7rv/t2 - 1/4) - l)dt = 
A/2 960 



we obtain in the same way as before 

" Sign(t)6(t) " °' 

where G(t) = F(t) - sign(t) I |^t 4 . 

Again the Fourier Tauberian theorem in |Saf01] implies 

Theorem E.4. The function satisfies 

1 , , 1 4 2Z> 2 + 7TZ> / D \ 3 11 

N l (r) < r 4 H r H H 

xK ' ~ 16tt 4tt 2 c2(x) V d(x)7 87T15 

1 „ i> 2 / z> \ 3 11 



Nl(r) > — r 4 - — — — r + 



16tt 27r 2 rf(x) V. J 8tt 15 



where ^ is the same as z/ 3 in |Saf01] and satisfies the esimtate z/ 3 < 6-^3 < 8. 
This results in the bound for the derivative of the eigenvalues 

Corollary E.5. Let ip(x) be an eigenfunction of the Laplace operator with eigenvalue 
A such that H^IU 2 = 1 and let n x a unit tangent vector at x G X . Then, 

, l2 4z> 2 + 7ri> / z> \ 3 11 
|m/>(x)| < 



4n 2 d(x) V d(x) J 4vr 15 
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Again, for numerical purposes we will simply use the bound 

£ fi ( r+ r) ,+ M- 
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